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ABSTRACT 

Accounting for the Reynolds number is critical in numerical simulations of turbulence, partic- 
ularly for subsonic flow. For Smoothed Particle Hydrodynamics (SPH) with constant artificial 
viscosity coefficient a, it is shown that the effective Reynolds number in the absence of ex- 
plicit physical viscosity terms scales linearly with the Mach number — compared to mesh 
schemes, where the effective Reynolds number is largely independent of the flow velocity. 
As a result, SPH simulations with a = 1 will have low Reynolds numbers in the subsonic 
regime compared to mesh codes, which may be insufficient to resolve turbulent flow. This 
explains the failure of |B auer & Springel] ( |2Q 1 1 1 arXivil 109.441 3vl) to find agreement be- 
tween the moving-mesh code AREPO and the gadget SPH code on simulations of driven, 
subsonic ~ O.Scg) turbulence appropriate to the intergalactic/intracluster medium, where it 
was alleged that SPH is somehow fundamentally incapable of producing a Kolmogorov-like 
turbulent cascade. We show that turbulent flow with a Kolmogorov spectrum can be easily 
recovered by employing standard methods for reducing a away from shocks. 
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^ 1 INTRODUCTION 

I Turbulence in astrophysics is of key importance for the interstellar 
(ISM), intracluster (ICM) and intergalactic medium (IGM). Com- 
pressible, hydrodynamic turbulence is characterised by two di- 
• 1^ mensionless parameters, the Mach number M = V/cs, and the 
^ Reynolds number fStokes|185l] |Reynolds|1883l ) 

c5 7^e^^, (1) 

where V is the flow velocity, L is a typical length scale, u is the 
viscosity of the fluid and Cs is the sound speed. Physically, these 
parameters estimate the relative importance of each of the terms 
in the Navier-Stokes equations — the Mach number specifies the 
ratio of the inertial term, (v • V)v, to the pressure term, VP/ p, 
while the Reynolds number specifies the ratio of the inertial term 
to the viscous dissipation term, z/V^v. Mathematically, these two 
parameters — along with the boundary conditions and driving — 
entirely characterise the flow. 

Given the importance of turbulence in theoretical models, it is 
crucial that agreement can be found between codes used for simula- 
tions of the ISM and ICM/IGM. Several comparison projects have 
been published recently comparing simulations of both decaying 
( |Kitsionas erari|2009l ) and driven { [Price & Federrath||2010a| l su- 
personic turbulence relevant to molecular clouds. However, fewer 
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calculations appropriate to the ICM or IGM have been performed. 
In a recent preprint [Bauer & Springel] ( |20 11 1 arXiv:1109.4413Yl) 
have set out to extend the high Mach number comparisons to the 
mildly compressible, driven, subsonic turbulence thought to be ap- 
propriate to the ICM and IGM. In this case the motions are com- 
parable to or smaller than the sound speed, turbulent motions are 
dissipated by viscosity, and the flow is mainly characterised by the 
Reynolds number, similar to turbulence in the laboratory. In par- 
ticular, it is well known from laboratory studies that the transition 
from laminar flow to fully developed turbulence only occurs once a 
critical Reynolds number is reached — for example for Poiseuille 
flow (water flowing in a pipe) this is observed for IZe > 2000 (e.g. 
[Reynolds] 1895] ). 

Since at low Mach number the Reynolds number controls not 
only the transition to turbulence but also the character of such tur- 
bulence (e.g. the extent of the inertial range) it is critical to specify, 
or at least estimate, the Reynolds number employed in numerical 
simulations of turbulence in order to compare with the physical 
Reynolds numbers in the problems of interest. For the ISM, the 
physical Reynolds numbers are high (e.g. jElmegreen & Scalo[2004| 
estimate IZe ^ 10^-10^ for the cold ISM) so the approach adopted 
has been to fix the Mach number and try to reach high numeri- 
cal Reynolds numbers by minimising numerical dissipation away 
from shocks. Estimates for IZe in the ICM/IGM are more difficult. 
[Brunetti & Lazarian] p007| ) estimate T^e ^ 52 but point out that 
extremely high values > 10^° are also possible in the presence of 
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magnetic fields, which change not only the viscosity but also intro- 
duce new scales into the problem (see |Lazarian|2006l ). This means 
that a range of studies are necessary, either by explicitly adding 
physical viscosity terms to reach a desired (low) IZe, or by seeking 
to minimise numerical viscosity to reach very high IZe. A minimum 
condition for production of turbulent flow in both simulations and 
reality is that the Reynolds numbers should at least be high enough 
for the development of turbulence (i.e. IZe ^ 10^). 

In their study, Bauer & SpringeT|pOlT| compare the moving- 
mesh code AREPO (run in both moving and fixed-grid modes) 
and the Smoothed Particle Hydrodynamics (SPH) code GADGET- 
3. The main conclusion of the paper is that, despite good results 
in the supersonic regime (in agreement with [Price & Federrath] 
|2Q10a| l, "the widely employed standard formulation of SPH fails 
quite badly in the subsonic regime", because of a failure to build 
up a "Kolmogorov-like turbulent cascade" which is apparent in the 
moving and fixed mesh calculations at the resolutions employed 
and is expected on theoretical grounds ( Kolmogorov 1941) to oc- 
cur at high Reynolds number. Bauer & Springel,(2011^ attribute the 
disagreement between the codes to "large errors in SPH's gradi- 
ent estimate and the associated subsonic velocity noise", producing 
"essentially unphysical results in the subsonic regime". This "casts 
doubt on the reliability of SPH for simulations of cosmic structure 
formation". However, the Reynolds number is neither fixed nor es- 
timated in any of the calculations. 

In this Letter show how the Reynolds numbers in SPH turbu- 
lence calculations can be determined (Sec. [2]). This suggests that 
in fact the main issue in the SPH calculations employed by [Bauer | 
|& Springe]] ( |2011| ) is that they employ constant — and thus large 
— viscosity parameters for their simulations, which in turn leads to 
low Reynolds numbers at subsonic velocities, explaining the failure 
to produce a Kolmogorov-like spectrum. To demonstrate this, we 
have rerun their cal culations (Sec. [3]), ado pting the standard SPH 
viscosity switch of [Morris & Monaghan| (T997| ) which increases 
the effective Reynolds number by roughly an order of magnitude 
and correspondingly leads to results in much better agreement with 
the grid-based simulations shown in their preprint and with Kol- 
mogorov's scaling relations. The results are discussed and sum- 
marised in Sec.|4] 



2 REYNOLDS NUMBERS IN SPH TURBULENCE 
CALCULATIONS 

An advantage of SPH, being a Hamiltonian method, is that the 
only dissipation present in the system is that which is explicitly 
added. This means that it is straightforward to derive Reynolds 
numbers for SPH calculations, since any terms added can be di- 
rectly translated into their physical equivalents. The standard SPH 
artificial viscosity term in 3D corresponds to Navier-Stokes viscos- 
ity terms with shear {v) and bulk (Q viscosity parameters given by 
(c.f. jLodato & Prrce]20l0l[Price|2012t 



(2) 



where h is the smoothing length and a is the SPH artificial viscos- 
ity parameter (Note that in 2D the parameters are 1/8 and 5 / 24, 
respectively, c.f. M urray] 1996[ ponaghan 2005 , Price,2012j . For 
low Mach number calculations we can expect that the V(V • v) 
terms represent only a small contribution to the overall dissipation 
rate, such that the dissipation is mainly controlled by the shear vis- 
cosity term. Furthermore the maximum signal velocity Vsig ~ Cs at 



low Mach number, such that the effective Reynolds number can be 
easily computed according to 



7^e = 



VL 



a n 



(3) 



where M is the Mach number. Note that, with A4 fixed, the 
Reynolds number is determined entirely by two parameters: the 
value of a and the numerical resolution h/L. For simulations in a 
periodic box using the parameters employed by [Bauer & Springe!] 
([2011j) we have 



7^e 2.4nx 
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0.3 
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64 



-1/3 



(4) 



where is the number of particles along the box length L and 
A^ngb is the neighbour number parameter in GADGET- 3, which cor- 
responds to 



iVngb,3D = -n{r]R)^, 



(5) 



where R is the kernel truncation radius in units of h (i.e. R — 2 
for the cubic spline kernel) and 77 is the usual parameter speci- 
fying the smoothing length in units of the particle spacing {h = 
?7[m/p]^/^^^™). Referring to the neighbour number is in general 
misleading since it is not independent of the resolution length h 
meaning that increasing A^ngb also corresponds to increasing h (see 
discussion in |Price|2012[ ). In this Letter we use rj — 1.2 correspond- 
ing to ^ 58 neighbours in a uniform particle distribution. 

From Eq.|4]it is easy to see why simulations are difficult at low 
Mach number {M < 0.5) in SPH with constant viscosity param- 
eters. For the resolutions employed by [Bauer & Spring'el|p011| ), 
namely rix — 64, 128 and 256 the Reynolds numbers are given 
by IZe = 154, IZe = 307 and IZe = 614 respectively. By com- 
parison, the Reynolds numbers that would be reached even with 
a = 1 at Mach 10 would be approximately 30 times highe r, giving 



IZe ^ 18, 000 at 2 56^ and 7^e ^ 37, 000 at 512^ particles. [Price & 
[Federrath (2010a) also employed viscosity switches that decrease 
the viscosity away from shocks by up to an order of magnitude, 
meaning that they obtain Reynolds numbers of ~ 10^ and higher 
in practice (a plot of the Reynolds number in the "Price & Federrath] 
|2010a calculations is shown in Fig. 2 of [Price & Federrath 2010b, 
Tie of up to 10^ are achieved in the densest regions where the SPH 
resolution is highest). 

Note also that the situation regarding numerical viscosity is 
very different in SPH compared to grid-based codes. The intrin- 
sic numerical dissipation in a grid-based code is in general lin- 
early proportional to the advection velocity relative to the grid (e.g. 
[Robertson et al.|20T0| , whereas in SPH the artificial viscosity term 
is linearly proportional to the resolution and largely independent of 
the velocity except at high Mach number where the /3 term comes 
into play. This means that Reynolds numbers are roughly constant 
in Eulerian schemes over a range of Mach numbers — because the 
increase in V is correspondingly offset by the increase in the nu- 
merical viscosity, whereas the Reynolds number in SPH has a linear 
Mach number dependence. Given that good agreement was found 
between grid and SPH codes in the power spectra in the [Price &[ 
[Federrath] \20 1 Oaj ) calculations, it is likely that the Reynolds num- 
ber achievable in the fixed and moving grid schemes employed in 
AREPO are of similar magnitude (i.e. ^ 10^ at 512^ and slightly 
lower at 256^). 
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Figure 1. Cross section slices through the midplane showing magnitud e of the velocity field (top ro w), density (second row) and |V x vp (bottom row) 
at f = 10 for SPH calculations with the phantom code employing the [Morris & Monaghan||T997) vi scosity switch at resolutio ns of 64^, 128^ and 256^ 
particles (left to right). The 256^ calculations may be compared to the corresponding panels in Fig. 3 of |Bauer & Springel||201l} . Estimates of the effective 
Reynolds numbers are indicated for each calculation. 



3 REDUCING THE VISCOSITY IN SPH 



3.1 Standard approaches 



[Bauer & S pringel (201 1 ) do not use any viscosity switches for their 
main calculations, despite the fact that most of these switches are 
at least ^15 years old and in widespread use. The standard vis- 
cosity switch in use is the one proposed by [Morris & MonaghaiT 
( |1997 |), where a is a time-dependent parameter that responds to a 
source term proportional to — V • v (i.e., converging flows) and in 
the absence of such terms decays to a minimum ttmin, typically set 
to 0.1. Already use of this switch would substantially increase the 
Reynolds number, though even a factor of 10 reduction in a gives 
only IZe ~ 6000 for their ''S3" calculation which is still a far cry 
from the 7^e ^ 10^ achieved by [Price & Federrath |( [2010a| . 



Fig. [T] shows the results of a series of 3 calculations performed 
with the PHANTOM SPH code ( [Price & Federrath[[2010a| |Lodato[ 
|& Price[[2010| ), employing the same driving routine as described 



in 'Price & Federrath[ (f2010a^ (adapted from'Federrath et al.'2008' 
[2010 j with the same par ameters as the Bauer & Springel ( 2011^ 
calculationQbut with the Morris & Monaghan ( 1997| ) switch with 
ttmin = 0.05, resulting in a ^ 0.1. The runs are performed 
in a periodic box x^y^z G [0, 1], using an isothermal equation 
of state with a sound speed in code units of unity. Although the 
Reynolds numbers achieved are evidently still lower than in the 



^ The SPH version of the driving routine can be made available on re- 
quest. The parameters used here are stenergy = 0.002, stdecay = 1-0, 

•S^solweight — 1.0, ststirmin — 6.28, sfgtirmax — 18.85, st^tfreq ~ 

0.005, corresponding to the stirring energy, decay timescale, solenoidal 
driving, minimum and maximum wavenumbers and the frequency with 
which the stirring is updated. The driving is given a k~^^^ wavenumber 
dependence in the stirring range as described in [Bauer & Springel] {2011} 
and denoted as stspectform = 2 in the input file. The random number gen- 
erator and seed will however differ from their calculations. 
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4 DISCUSSION 



Figure 2. Time-averaged fc^/^ -compensated power spectra from subsonic 
SPH turbulence calculations using the [Morris & MonaglianH1997'| viscos- 
ity switch at a resolution of 64^, 128^ and 256^ particles, as indicated, 
for which the corresponding Reynolds numbers are ~ 1500, 3000 and 
6000, respectively. The shaded regions show the Icr errors from the time- 
averaging. At the highest Reynolds numbers a Kolmogorov-like k~^^^ 
slope is evident at large scales. 



grid-based calculations employed by I Bauer & Springel] ( |2011| ), it 
is clear that already this is a dramatic improvement on the SPH 
simulations shown in their preprint. The resulting power spectra 
are shown in Fig. [2] showing the time-averaged spectrum from 191 
snapshots sampled every At = 0.1 between t = 6 and t = 25. The 
^— axis shows the power spectrum compensated by k^^^ such that 
a k~^^^ spectrum would appear horizontal. Though the spectrum 
"turns over" at relatively small k at low resolution, a clear k~^^^ 
range is apparent in the highest resolution calculations. The resolu- 
tion dependence of the high k turnover is also consistent with the 
expected IZe^^^ dependence of the dissipation scale. 



3.2 The state of the art 

With mean motions that are around 1/3 of the sound speed and 
transsonic fluctuations, one cannot simply reduce the SPH artificial 
viscosity parameters arbitrarily (such as the a = 0.1 and a = 
calculations attempted by [Bauer & Springel|2011| ), since this term 
is necessary to capture the physical dissipation that occurs due to 
the non-linear steepening of waves. Such an approach may be ade- 
quate for very low Mach number (i.e., incompressible) calculations 
but it provides no easy answer at ^ 0.3cs. Instead it is clear that 
to achieve similar results with SPH an improved viscosity switch is 
necessary in order to both capture non-linear steepening and shocks 
as well as reducing the viscosity to very low values where it is not 
needed. The switch proposed recently by jCuUen & Dehnen| ( |2010| ) 
represents the current state-of-the-art in this regard, essentially a 
thoroughly enhanced and improved version of the [Morris & MorT] 
|aghan|(r997| ) approach. In particular, they show that they are able 
to simulate linear waves for over 50 periods with essentially no nu- 
merical dissipation, using the same parameters as would be applied 
in shock problems. Thus, with an implementation of the jCuUen &| 
|Dehnen| ( |2010| ) switch it may be expected that significantly higher 
Reynolds numbers are achievable in SPH. 



[Bauer & Springel| p011| ) argue that "large errors in SPH's gra- 
dient estimate" are responsible for the failure to reproduce a 
Kolmogorov-like turbulent cascade in their SPH calculations. Fig- 
ure [2| demonstrates that this argument is incorrect, since we are 
able to obtain a k~^^^ spectrum using only standard SPH gradi- 
ent terms and a very similar SPH neighbour number to that em- 
ployed in their preprint. However, we find that the appearance of 
a power-law inertial range in the power spectrum strongly depends 
on the Reynolds number employed in the calculations, requiring 
at least IZe ^ 1500. This explains the failure to produce a tur- 
bulent cascade in their SPH results, since the maximum Reynolds 
numbers they achieve are ^ 600. With the [Morris & Monaghanj 
( |1997f viscosity switch employed in this Letter we estimate that 
we are able to achieve T^e ~ 6, 000 at 256^ particles which al- 
ready brings the SPH results into much better agreement with the 
grid-based results shown in their work. Indeed, both [Dolag et"aL] 
(2005 ) and Vald arninilp011| ) have already pointed out that using 
this switch could substantially improve SPH simulations of turbu- 
lence in galaxy clusters. 

It should be noted that |Bauer & Springel|p01 1) do experiment 
with reduced viscosity parameters in their preprint, using either a 
fixed a = 0.1 or the Balsara (1995 ) switch and also a run with 
zero viscosity (as we have already discussed in Sec. [3] it is not 
clear that one can simply reduce the parameters arbitrarily, so this 
approach is questionable — particularly the a = calculation). 
Indeed, both the a = 0.1 and Balsara- switch calculations show a 
dramatic improvement in the power spectrum at large scales. The 
authors dismiss this result because of a corresponding increase in 
power at /c > 100. However, the power at these scales is low ampli- 
tude {r^ 10~^) and thus sensitive to all manner of numerical arte- 
facts (e.g. the interpolation procedure as demonstrated in Fig. 4 of 
their preprint). Indeed, we do not find the upturn in power at large 
k seen in their results (c.f. Fig. [2]), most likely due to our improved 
power spectrum estimation — here computed by interpolating the 
SPH data to a 3D grid using the kernel and employing a Fast Fourier 
Transform, rather than the "nearest neighbour sampling" procedure 
employed in their preprint. 

Finally, it is important to compare the Reynolds numbers 
achievable in numerical simulations to the physical Reynolds num- 
bers in the problems of interest (c.f. Sec. [TJ. For the cold ISM, 
IZe ^ 10^-10^ which, though high, is not as high as is often as- 
sumed, and is certainly within reach of being resolved with cur- 
rently achievable resolutions (c.f. Sec. [2]). This implies at the very 
least that physical viscosity should be introduced into ISM turbu- 
lence simulations in the near future. In the ICM/IGM, IZe ~ 52 
would seem to imply laminar flow, though very high estimates for 
'^e (> 10^°) apply in the presence of magnetic fields. Reaching 
such Reynolds numbers is not presently achievable with any nu- 
merical code. However, this may also imply that ultimately it is 
quite incorrect to try to simulate purely hydrodynamic ICM/IGM 
turbulence at high Reynolds number without taking into account 
more detailed physics, such as magnetic fields. 



5 CONCLUSIONS 

In this Letter we have emphasized the importance of accounting for 
the Reynolds number in numerical turbulence simulations, partic- 
ularly in the subsonic regime where it is the main parameter con- 
trolling not only whether the flow is turbulent but also the char- 
acter of such turbulence. In particular differences in the Reynolds 



number can produce strong physical differences between calcula- 
tions, independent of any numerical influences. We have shown that 
use of viscosity switches is the key to higher Reynolds numbers 
in SPH turbulence calculations, producing spectra consistent with 
Kolmogorov theory. 
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